Right Censoring

• Patients are in a clinical trial for cancer, some on a new treatment and some on standard of care.
• Some patients in each group have died by the end of the study. We know the survival time (measured for example from time of diagnosis—each person on their own clock).
• Patients still alive at the end of the study are right censored.
• Patients who are lost to follow-up or withdraw from the study may be right-censored.

Left and Interval Censoring

• An individual tests positive for HIV.
• If the event is infection with HIV, then we only know that it has occurred before the testing time \(t\), so this is left censored.
• If an individual has a negative HIV test at time \(t_1\) and a positive HIV test at time \(t_2\), then the infection event is interval censored.

The Cumulative Distribution Function (CDF)

The cumulative distribution function is defined as:

\[ \begin{aligned} F(t) &\stackrel{\text{def}}{=}\Pr(T \le t)\\ &=\int_{u=-\infty}^t f(u) du \end{aligned} \]

Here are some examples of exponential cdfs:

The Survival Function

We can also view the hazard function as the derivative of the negative of the logarithm of the survival function:

The Cumulative Hazard Function

Since \({\lambda}(t) = \frac{\partial}{\partial t}{\left\{-\log{\text{S}(t)}\right\}}\) (see Theorem 6), we also have:

As we will see below, \({\Lambda}(t)\) is tractable to estimate, and we can then derive an estimate of the hazard function using an approximate derivative of the estimated cumulative hazard.

Some Key Mathematical Relationships among Survival Concepts

Diagram:

\[ \text{f}(t) \xleftarrow[\text{S}(t){\lambda}(t)]{-S'(t)} \text{S}(t) \xleftarrow[]{\text{exp}{\left\{-{\Lambda}(t)\right\}}} {\Lambda}(t) \xleftarrow[]{\int_{u=0}^t {\lambda}(u)du} {\lambda}(t) \xleftarrow[]{\text{exp}{\left\{\eta(t)\right\}}} \eta(t) \]

\[ \text{f}(t) \xrightarrow[\int_{u=t}^\infty \text{f}(u)du]{\text{f}(t)/{\lambda}(t)} \text{S}(t) \xrightarrow[-\log{\text{S}(t)}]{} {\Lambda}(t) \xrightarrow[{\Lambda}'(t)]{} {\lambda}(t) \xrightarrow[\text{log}{\left\{{\lambda}(t)\right\}}]{} \eta(t) \]

Identities:

\[ \begin{aligned} \text{S}(t) &= 1 - \text{F}(t)\\ &= \text{exp}{\left\{-{\Lambda}(t)\right\}}\\ \text{S}'(t) &= -f(t)\\ {\Lambda}(t) &= -\text{log}{\left\{\text{S}(t)\right\}}\\ {\Lambda}'(t) &= {\lambda}(t)\\ {\lambda}(t) &= \frac{\text{f}(t)}{\text{S}(t)}\\ &= -\frac{\partial}{\partial t}\log{\text{S}(t)} \\ f(t) &= {\lambda}(t)\cdot \text{S}(t)\\ \end{aligned} \]

Some proofs (others left as exercises):

\[ \begin{aligned} S'(t) &= \frac{\partial}{\partial t}(1-F(t))\\ &= -F'(t)\\ &= -f(t)\\ \end{aligned} \]

\[ \begin{aligned} \frac{\partial}{\partial t}\log{\text{S}(t)} &= \frac{S'(t)}{\text{S}(t)}\\ &= -\frac{f(t)}{\text{S}(t)}\\ &= -{\lambda}(t)\\ \end{aligned} \]

\[ \begin{aligned} {\Lambda}(t) &\stackrel{\text{def}}{=}\int_{u=0}^t h(u) du\\ &= \int_0^t -\frac{\partial}{\partial u}\text{log}\left\{S(u)\right\} du\\ &= \left[-\text{log}\left\{S(u)\right\}\right]_{u=0}^{u=t}\\ &= \left[\text{log}\left\{S(u)\right\}\right]_{u=t}^{u=0}\\ &= \text{log}{\left\{\text{S}(0)\right\}} - \text{log}{\left\{\text{S}(t)\right\}}\\ &= \text{log}{\left\{1\right\}} - \text{log}{\left\{\text{S}(t)\right\}}\\ &= 0 - \text{log}{\left\{\text{S}(t)\right\}}\\ &=-\text{log}{\left\{\text{S}(t)\right\}} \end{aligned} \]

Corollary:

\[\text{S}(t) = \text{exp}\left\{-{\Lambda}(t)\right\}\]

Example: Time to death the US in 2004

The first day is the most dangerous:

[R code]

# download `survexp.rda` from:
# paste0(
# "https://github.com/therneau/survival/raw/",
# "f3ac93704949ff26e07720b56f2b18ffa8066470/",
# "Data/survexp.rda")

# (newer versions of `survival` don't have the first-year breakdown; see:
# https://cran.r-project.org/web/packages/survival/news.html)

fs::path(
  here::here(),
  "Data",
  "survexp.rda"
) |>
  load()
s1 <- survexp.us[, "female", "2004"]
age1 <- c(
  0.5 / 365.25,
  4 / 365.25,
  17.5 / 365.25,
  196.6 / 365.25,
  1:109 + 0.5
)
s2 <- 365.25 * s1[5:113]
s2 <- c(s1[1], 6 * s1[2], 22 * s1[3], 337.25 * s1[4], s2)
cols <- rep(1, 113)
cols[1] <- 2
cols[2] <- 3
cols[3] <- 4

plot(age1, s1, type = "b", lwd = 2, xlab = "Age", ylab = "Daily Hazard Rate", col = cols)

text(10, .003, "First Day", col = 2)
text(18, .00030, "Rest of First Week", col = 3)
text(18, .00015, "Rest of First month", col = 4)

Likelihood with censoring

If an event time \(T\) is observed exactly as \(T=t\), then the likelihood of that observation is just its probability density function:

\[ \begin{aligned} \mathscr{L}(t) &= \text{f}(T=t)\\ &\stackrel{\text{def}}{=}\text{f}_T(t)\\ &= {\lambda}_T(t)\text{S}_T(t)\\ \ell(t) &\stackrel{\text{def}}{=}\text{log}{\left\{\mathscr{L}(t)\right\}}\\ &= \text{log}{\left\{{\lambda}_T(t) \text{S}_T(t)\right\}}\\ &= \text{log}{\left\{{\lambda}_T(t)\right\}} + \text{log}{\left\{\text{S}_T(t)\right\}}\\ &= \text{log}{\left\{{\lambda}_T(t)\right\}} - {\Lambda}_T(t)\\ \end{aligned} \]

If instead the event time \(T\) is censored and only known to be after time \(y\), then the likelihood of that censored observation is instead the survival function evaluated at the censoring time:

\[ \begin{aligned} \mathscr{L}(y) &=p_T(T>y)\\ &\stackrel{\text{def}}{=}\text{S}_T(y)\\ \ell(y) &\stackrel{\text{def}}{=}\text{log}{\left\{\mathscr{L}(y)\right\}}\\ &=\text{log}{\left\{\text{S}(y)\right\}}\\ &=-{\Lambda}(y)\\ \end{aligned} \]

\[ \begin{aligned} \mathscr{L}(y,d) &= \text{p}(Y=y, D=d)\\ &= {\left[\text{p}(T=y,C> y)\right]}^d \cdot {\left[\text{p}(T>y,C=y)\right]}^{1-d}\\ \end{aligned} \]

\[ \begin{aligned} \mathscr{L}(y,d) &= {\left[\text{p}(T=y,C> y)\right]}^d\cdot {\left[\text{p}(T>y,C=y)\right]}^{1-d}\\ &= {\left[\text{p}(T=y)\text{p}(C> y)\right]}^d\cdot {\left[\text{p}(T>y)\text{p}(C=y)\right]}^{1-d}\\ &= {\left[\text{f}_T(y)\text{S}_C(y)\right]}^d\cdot {\left[\text{S}(y)\text{f}_C(y)\right]}^{1-d}\\ &= {\left[\text{f}_T(y)^d \text{S}_C(y)^d\right]}\cdot {\left[\text{S}_T(y)^{1-d} \text{f}_C(y)^{1-d}\right]}\\ &= {\left(\text{f}_T(y)^d \cdot \text{S}_T(y)^{1-d}\right)}\cdot {\left(\text{f}_C(y)^{1-d} \cdot \text{S}_C(y)^{d}\right)} \end{aligned} \]

\[ \begin{aligned} \ell(y,d) &= \text{log}{\left\{\mathscr{L}(y,d)\right\}}\\ &= \text{log}{\left\{ {\left(f_T(y)^d \cdot S_T(y)^{1-d}\right)} \cdot {\left(f_C(y)^{1-d} \cdot S_C(y)^{d}\right)} \right\}} \\ &= \text{log}{\left\{f_T(y)^d \cdot S_T(y)^{1-d}\right\}} + \text{log}{\left\{f_C(y)^{1-d} \cdot S_C(y)^{d}\right\}} \end{aligned} \] Let

• \(\theta_T\) represent the parameters of \(p_T(t)\),
• \(\theta_C\) represent the parameters of \(p_C(c)\),
• \(\theta = (\theta_T, \theta_C)\) be the combined vector of all parameters.

\[ \begin{aligned} \ell'(y,d) &= \frac{\partial}{\partial \theta} {\left[ \text{log}{\left\{f_T(y)^d \cdot S_T(y)^{1-d}\right\}} + \text{log}{\left\{f_C(y)^{1-d} \cdot S_C(y)^{d}\right\}} \right]} \\ &= {\left( \frac{\partial}{\partial \theta} \text{log}{\left\{ f_T(y)^d \cdot S_T(y)^{1-d} \right\}} \right)} + {\left( \frac{\partial}{\partial \theta} \text{log}{\left\{ f_C(y)^{1-d} \cdot S_C(y)^{d} \right\}} \right)} \end{aligned} \]

\[ \begin{aligned} \ell'_{\theta_T}(y,d) &\stackrel{\text{def}}{=}\frac{\partial}{\partial \theta_T}\ell(y,d)\\ &= \left( \frac{\partial}{\partial \theta_T} \text{log}\left\{ f_T(y)^d \cdot S_T(y)^{1-d} \right\} \right) + \left( \frac{\partial}{\partial \theta_T} \text{log}\left\{ f_C(y)^{1-d} \cdot S_C(y)^{d} \right\} \right)\\ &= \left( \frac{\partial}{\partial \theta_T} \text{log}\left\{ f_T(y)^d \cdot S_T(y)^{1-d} \right\} \right) + 0\\ &= \frac{\partial}{\partial \theta_T} \text{log}\left\{ f_T(y)^d \cdot S_T(y)^{1-d} \right\}\\ \end{aligned} \]

Then:

\[ \begin{aligned} \mathscr{L}(y,d) &= f_T(y)^d \cdot S_T(y)^{1-d}\\ &= \left(h_T(y)^d S_T(y)^d\right) \cdot S_T(y)^{1-d}\\ &= h_T(y)^d \cdot S_T(y)^d \cdot S_T(y)^{1-d}\\ &= h_T(y)^d \cdot S_T(y)\\ &= S_T(y) \cdot h_T(y)^d \\ \end{aligned} \]

\[ \begin{aligned} \ell(y,d) &=\text{log}\left\{\mathscr{L}(y,d)\right\}\\ &= \text{log}\left\{S_T(y) \cdot h_T(y)^d\right\}\\ &= \text{log}\left\{S_T(y)\right\} + \text{log}\left\{h_T(y)^d\right\}\\ &= \text{log}\left\{S_T(y)\right\} + d\cdot \text{log}\left\{h_T(y)\right\}\\ &= -H_T(y) + d\cdot \text{log}\left\{h_T(y)\right\}\\ \end{aligned} \]

Exponential Distribution

• The exponential distribution is the base distribution for survival analysis.
• The distribution has a constant hazard \(\lambda\)
• The mean survival time is \(\lambda^{-1}\)

Mathematical details of exponential distribution

\[ \begin{aligned} f(t) &= \lambda \text{e}^{-\lambda t} \\ F(t) &= 1-\text{e}^{-\lambda t}\\ \text{S}(t)&= \text{e}^{-\lambda t}\\ \ln(\text{S}(t))&=-\lambda t\\ {\lambda}(t) &= -\frac{f(t)}{\text{S}(t)} = -\frac{\lambda \text{e}^{-\lambda t}}{\text{e}^{-\lambda t}}=\lambda \\ E(t) &= \lambda^{-1} \\ Var(t) &= \lambda^{-2} \\ \text{log}{\left\{f(t)\right\}} &= \text{log}{\left\{\lambda\right\}} - \lambda t \\ \frac{\partial}{\partial \lambda} \text{log}{\left\{f(t)\right\}} &= \lambda^{-1} - t \\ &= \text{E}{\left[t\right]} - t \\ &= -(\text{E}{\left[t\right]} - t) \\ &= -\varepsilon \end{aligned} \]

Prediction intervals for time-to-event outcomes

Estimating \(\lambda\)

• Suppose we have \(m\) exponential survival times of \(t_1, t_2,\ldots,t_m\) and \(k\) right-censored values at \(u_1,u_2,\ldots,u_k\).

• A survival time of \(t_i=10\) means that subject \(i\) died at time 10. A right-censored time \(u_i=10\) means that at time 10, subject \(i\) was still alive and that we have no further follow-up.

• For the moment we will assume that the survival distribution is exponential and that all the subjects have the same parameter \(\lambda\).

We have \(m\) exponential survival times of \(t_1, t_2,\ldots,t_m\) and \(k\) right-censored values at \(u_1,u_2,\ldots,u_k\). The log-likelihood of an observed survival time \(t_i\) is \[ \text{log}\left\{\lambda \text{e}^{-\lambda t_i}\right\} = \text{log}\left\{\lambda\right\}-\lambda t_i \] and the likelihood of a censored value is the probability of that outcome (survival greater than \(u_j\)) so the log-likelihood is

\[ \begin{aligned} \ell_j(\lambda) &= \text{log}\left\{\lambda \text{e}^{u_j}\right\} \\ &= -\lambda u_j \end{aligned} \]

\[ \begin{aligned} \ell''&=-m/\lambda^2\\ &< 0\\ \hat E[T] &= \hat\lambda^{-1}\\ &= \frac{T+U}{m} \end{aligned} \]

Fisher Information and Standard Error

\[ \begin{aligned} E[-\ell''] & = m/\lambda^2\\ \text{Var}{\left(\hat\lambda\right)} &\approx {\left(E[-\ell'']\right)}^{-1}\\ &=\lambda^2/m\\ \text{SE}\left(\hat\lambda\right) &= \sqrt{\text{Var}{\left(\hat\lambda\right)}}\\ &\approx \lambda/\sqrt{m} \end{aligned} \]

Other Parametric Survival Distributions

• Any density on \([0,\infty)\) can be a survival distribution, but the most useful ones are all skew right.
• The most frequently used generalization of the exponential is the Weibull.
• Other common choices are the gamma, log-normal, log-logistic, Gompertz, inverse Gaussian, and Pareto.
• Most of what we do going forward is non-parametric or semi-parametric, but sometimes these parametric distributions provide a useful approach.

Basic ideas

• Mostly, we work without a parametric model.

• The first task is to estimate a survival function from data listing survival times, and censoring times for censored data.

• For example one patient may have relapsed at 10 months. Another might have been followed for 32 months without a relapse having occurred (censored).

• The minimum information we need for each patient is a time and a censoring variable which is 1 if the event occurred at the indicated time and 0 if this is a censoring time.

Data documentation for drug6mp

# library(printr) # inserts help-file output into markdown output
library(KMsurv)
?drug6mp

Exponential model

\[\hat\lambda = \frac{\sum_{i=1}^nD_i}{\sum_{i=1}^nY_i}\]

\[ \begin{aligned} \hat \lambda_{\text{placebo}} &= \frac{\sum_{i=1}^nD_i}{\sum_{i=1}^nY_i} \\ &= \frac{\sum_{i=1}^n1}{\sum_{i=1}^nY_i} \\ &= \frac{n}{\sum_{i=1}^nY_i} \\ &= \frac{1}{\bar{Y}} \\ &= \frac{1}{8.666667} \\ &= 0.115385 \end{aligned} \]

• For the 6-MP group, \(\hat\lambda = 9/359 = 0.025\)

\[ \begin{aligned} \hat \lambda_{\text{6-MP}} &= \frac{\sum_{i=1}^nD_i}{\sum_{i=1}^nY_i} \\ &= \frac{9}{359} \\ &= 0.02507 \end{aligned} \]

• The estimated hazard in the placebo group is 4.6 times as large as in the 6-MP group (assuming the hazard is constant over time).

Estimating survival in datasets without censoring

Estimating survival in datasets with censoring

• For the 6-MP patients, we cannot ignore the censored data because we know that the time to relapse is greater than the censoring time.

• For any time \(t\) in months, we know that 6-MP patients with times greater than \(t\) have not relapsed, and those with relapse time less than \(t\) have relapsed, but we don’t know if patients with censored time less than \(t\) have relapsed or not.

• The procedure we usually use is the Kaplan-Meier product-limit estimator of the survival function.

• The Kaplan-Meier estimator is a step function (like the empirical cdf), which changes value only at the event times, not at the censoring times.

• At each event time \(t\), we compute the at-risk group size \(Y\), which is all those observations whose event time or censoring time is at least \(t\).

• If \(d\) of the observations have an event time (not a censoring time) of \(t\), then the group of survivors immediately following time \(t\) is reduced by the fraction \[\frac{Y-d}{Y}=1-\frac{d}{Y}\]

Kaplan-Meier curve for drug6mp data

Here is the Kaplan-Meier estimated survival curve for the patients who received 6-MP in the drug6mp dataset (we will see code to produce figures like this one shortly):

[R code]

# | echo: false

require(KMsurv)
data(drug6mp)
library(dplyr)
library(survival)

drug6mp_km_model1 <-
  drug6mp |>
  mutate(surv = Surv(t2, relapse)) |>
  survfit(formula = surv ~ 1, data = _)

library(ggfortify)
drug6mp_km_model1 |>
  autoplot(
    mark.time = TRUE,
    conf.int = FALSE
  ) +
  expand_limits(y = 0) +
  xlab("Time since diagnosis (months)") +
  ylab("KM Survival Curve")

Summary

For the 6-MP patients at time 6 months, there are 21 patients at risk. At \(t=6\) there are 3 relapses and 1 censored observations.

The Kaplan-Meier factor is \((21-3)/21 = 0.857\). The number at risk for the next time (\(t=7\)) is \(21-3-1=17\).

At time 7 months, there are 17 patients at risk. At \(t=7\) there is 1 relapse and 0 censored observations. The Kaplan-Meier factor is \((17-1)/17 = 0.941\). The Kaplan Meier estimate is \(0.857\times0.941=0.807\). The number at risk for the next time (\(t=9\)) is \(17-1=16\).

Now, let’s graph this estimated survival curve using ggplot():

[R code]

library(ggplot2)
conflicts_prefer(dplyr::filter)
km.6mp |>
  ggplot(aes(x = t2, y = `KM Survival Curve`)) +
  geom_step() +
  geom_point(data = km.6mp |> filter(Censored > 0), shape = 3) +
  expand_limits(y = c(0, 1), x = 0) +
  xlab("Time since diagnosis (months)") +
  ylab("KM Survival Curve") +
  scale_y_continuous(labels = scales::percent)

The Surv function

To use the survival package, the first step is telling R how to combine the exit time and exit reason (censoring versus event) columns. The Surv() function accomplishes this task.

Example: Surv() with drug6mp data

library(survival)
drug6mp.v3 <-
  drug6mp.v2 |>
  mutate(
    surv2 = Surv(
      time = t2,
      event = (outcome == "relapsed")
    )
  )

print(drug6mp.v3)
#> # A tibble: 21 × 7
#>     pair remstat     t1    t2 relapse outcome   surv2
#>    <int> <chr>    <int> <int>   <int> <chr>    <Surv>
#>  1     1 partial      1    10       1 relapsed    10 
#>  2     2 complete    22     7       1 relapsed     7 
#>  3     3 complete     3    32       0 censored    32+
#>  4     4 complete    12    23       1 relapsed    23 
#>  5     5 complete     8    22       1 relapsed    22 
#>  6     6 partial     17     6       1 relapsed     6 
#>  7     7 complete     2    16       1 relapsed    16 
#>  8     8 complete    11    34       0 censored    34+
#>  9     9 complete     8    32       0 censored    32+
#> 10    10 complete    12    25       0 censored    25+
#> # ℹ 11 more rows

The output of Surv() is a vector of objects with class Surv. When we print this vector:

• observations where the event was observed are printed as the event time (for example, surv2 = 10 on line 1)

• observations where the event was right-censored are printed as the censoring time with a plus sign (+; for example, surv2 = 32+ on line 3).

The survfit function

Once we have constructed our Surv variable, we can calculate the Kaplan-Meier estimate of the survival curve using the survfit() function.

Example: survfit() with drug6mp data

Here we use survfit() to create a survfit object, which contains the Kaplan-Meier estimate:

drug6mp.km_model <- survfit(
  formula = surv2 ~ 1,
  data = drug6mp.v3
)

print.survfit() just gives some summary statistics:

print(drug6mp.km_model)
#> Call: survfit(formula = surv2 ~ 1, data = drug6mp.v3)
#> 
#>       n events median 0.95LCL 0.95UCL
#> [1,] 21      9     23      16      NA

summary.survfit() shows us the underlying Kaplan-Meier table:

summary(drug6mp.km_model)
#> Call: survfit(formula = surv2 ~ 1, data = drug6mp.v3)
#> 
#>  time n.risk n.event survival std.err lower 95% CI upper 95% CI
#>     6     21       3    0.857  0.0764        0.720        1.000
#>     7     17       1    0.807  0.0869        0.653        0.996
#>    10     15       1    0.753  0.0963        0.586        0.968
#>    13     12       1    0.690  0.1068        0.510        0.935
#>    16     11       1    0.627  0.1141        0.439        0.896
#>    22      7       1    0.538  0.1282        0.337        0.858
#>    23      6       1    0.448  0.1346        0.249        0.807

We can specify which time points we want using the times argument:

summary(
  drug6mp.km_model,
  times = c(0, drug6mp.v3$t2)
)
#> Call: survfit(formula = surv2 ~ 1, data = drug6mp.v3)
#> 
#>  time n.risk n.event survival std.err lower 95% CI upper 95% CI
#>     0     21       0    1.000  0.0000        1.000        1.000
#>    10     15       1    0.753  0.0963        0.586        0.968
#>     7     17       1    0.807  0.0869        0.653        0.996
#>    32      4       0    0.448  0.1346        0.249        0.807
#>    23      6       1    0.448  0.1346        0.249        0.807
#>    22      7       1    0.538  0.1282        0.337        0.858
#>     6     21       3    0.857  0.0764        0.720        1.000
#>    16     11       1    0.627  0.1141        0.439        0.896
#>    34      2       0    0.448  0.1346        0.249        0.807
#>    32      4       0    0.448  0.1346        0.249        0.807
#>    25      5       0    0.448  0.1346        0.249        0.807
#>    11     13       0    0.753  0.0963        0.586        0.968
#>    20      8       0    0.627  0.1141        0.439        0.896
#>    19      9       0    0.627  0.1141        0.439        0.896
#>     6     21       3    0.857  0.0764        0.720        1.000
#>    17     10       0    0.627  0.1141        0.439        0.896
#>    35      1       0    0.448  0.1346        0.249        0.807
#>     6     21       3    0.857  0.0764        0.720        1.000
#>    13     12       1    0.690  0.1068        0.510        0.935
#>     9     16       0    0.807  0.0869        0.653        0.996
#>     6     21       3    0.857  0.0764        0.720        1.000
#>    10     15       1    0.753  0.0963        0.586        0.968

Plotting estimated survival functions

We can plot survfit objects with plot(), autoplot(), or ggsurvplot():

library(ggfortify)
autoplot(drug6mp.km_model)

# not shown:
# plot(drug6mp.km_model)

# library(survminer)
# ggsurvplot(drug6mp.km_model)

quantiles of survival curve

We can extract quantiles with quantile():

drug6mp.km_model |>
  quantile(p = c(.25, .5)) |>
  as_tibble() |>
  mutate(p = c(.25, .5)) |>
  relocate(p, .before = everything())

For \(t \in t_1,...,t_n\): \[ \begin{aligned} \hat{{\lambda}}_t &= \frac{\sum_x m_{x,t}}{\sum_x n_{x,t}} \\ \hat{E}_{t,x} &= \hat{{\lambda}}_t * n_{x,t} \end{aligned} \]

The survdiff function

?survdiff

Example: survdiff() with drug6mp data

Now we are going to compare the placebo and 6-MP data. We need to reshape the data to make it usable with the standard survival workflow:

library(survival)
library(tidyr)
drug6mp.v4 <-
  drug6mp.v3 |>
  select(pair, remstat, t1, t2, outcome) |>
  # here we are going to change the data from a wide format to long:
  pivot_longer(
    cols = c(t1, t2),
    names_to = "treatment",
    values_to = "exit_time"
  ) |>
  mutate(
    treatment = treatment |>
      case_match(
        "t1" ~ "placebo",
        "t2" ~ "6-MP"
      ),
    outcome = if_else(
      treatment == "placebo",
      "relapsed",
      outcome
    ),
    surv = Surv(
      time = exit_time,
      event = (outcome == "relapsed")
    )
  )

Using this long data format, we can fit a Kaplan-Meier curve for each treatment group simultaneously:

drug6mp.km_model2 <-
  survfit(
    formula = surv ~ treatment,
    data = drug6mp.v4
  )

We can plot the curves in the same graph:

drug6mp.km_model2 |> autoplot()

We can also perform something like a t-test, where the null hypothesis is that the curves are the same:

o_e <- drug6mp.v4 |>
  arrange(exit_time) |>
  mutate(
    .by = treatment,
    n_exited = row_number(),
    n_at_risk = n() - n_exited + 1
  ) |>
  dplyr::summarize(
    .by = all_of(c("exit_time", "treatment")),
    n_at_risk = max(n_at_risk),
    n_events = sum(outcome == "relapsed")
  ) |>
  tidyr::pivot_wider(
    names_from = "treatment",
    values_from = c(n_at_risk, n_events)
  ) |>
  tidyr::fill(
    starts_with("n_at_risk"),
    .direction = "up"
  ) |>
  replace_na(list("n_events_placebo" = 0,
                  "n_events_6-MP" = 0)) |>
  mutate(
    n_at_risk = rowSums(across(starts_with("n_at_risk"))),
    n_events = rowSums(across(starts_with("n_events"))),
    marginal_hazard = n_events / n_at_risk,
    expected_6mp = marginal_hazard * `n_at_risk_6-MP`,
    expected_plc = marginal_hazard * n_at_risk_placebo,
    diff_6mp = `n_events_6-MP` - expected_6mp,
    diff_plc = n_events_placebo - expected_plc
  ) |>
  filter(n_events > 0)

o_e

o_e_summ <- o_e |>
  summarize(
    across(starts_with("expected"), sum),
    across(starts_with("n_events_"), sum)
  )
pander::pander(o_e_summ)

\[ \text{Var}{\left(O_i - E_i\right)} = \sum_j \frac {n_{1j} n_{2j} (m_j) (n_j - m_j)} {(n_j)^2 (n_j - 1)} \]

\[X^2 \approx \sum_{i = 1}^p \frac{(O_i - E_i)^2}{E_i}\]

with(
  o_e_summ,
  tibble(
    "6mp" = (`n_events_6-MP` - expected_6mp)^2 / expected_6mp,
    "placebo" = (n_events_placebo - expected_plc)^2 / expected_plc,
    sum = `6mp` + placebo
  )
) |>
  pander::pander()

R gives us both the exact and approximate results:

survdiff(
  formula = surv ~ treatment,
  data = drug6mp.v4
)
#> Call:
#> survdiff(formula = surv ~ treatment, data = drug6mp.v4)
#> 
#>                    N Observed Expected (O-E)^2/E (O-E)^2/V
#> treatment=6-MP    21        9     19.3      5.46      16.8
#> treatment=placebo 21       21     10.7      9.77      16.8
#> 
#>  Chisq= 16.8  on 1 degrees of freedom, p= 4e-05

Data from Copelan et al. (1991)

Study design

Treatment

• allogeneic (from a donor) bone marrow transplant therapy

Inclusion criteria

• acute myeloid leukemia (AML)
• acute lymphoblastic leukemia (ALL).

Possible intermediate events

• graft vs. host disease (GVHD): an immunological rejection response to the transplant
• platelet recovery: a return of platelet count to normal levels.

One or the other, both in either order, or neither may occur.

End point events

• relapse of the disease
• death

Any or all of these events may be censored.

KMsurv::bmt data in R

library(KMsurv)
?bmt

Analysis plan

• We concentrate for now on disease-free survival (t2 and d3) for the three risk groups, ALL, AML Low Risk, and AML High Risk.
• We will construct the Kaplan-Meier survival curves, compare them, and test for differences.
• We will construct the cumulative hazard curves and compare them.
• We will estimate the hazard functions, interpret, and compare them.

Survival Function Estimate and Variance

\[\hat{\text{S}}(t) = \prod_{t_i < t}\left[1-\frac{d_i}{Y_i}\right]\] where \(Y_i\) is the group at risk at time \(t_i\).

The estimated variance of \(\hat{\text{S}}(t)\) is:

We can use Equation 3 for confidence intervals for a survival function or a difference of survival functions.

Kaplan-Meier survival curves

library(KMsurv)
library(survival)
data(bmt)

bmt <-
  bmt |>
  as_tibble() |>
  mutate(
    group =
      group |>
      factor(
        labels = c("ALL", "Low Risk AML", "High Risk AML")
      ),
    surv = Surv(t2, d3)
  )

km_model1 <- survfit(
  formula = surv ~ group,
  data = bmt
)

library(ggfortify)
autoplot(
  km_model1,
  conf.int = TRUE,
  ylab = "Pr(disease-free survival)",
  xlab = "Time since transplant (days)"
) +
  theme_bw() +
  theme(legend.position = "bottom")

Understanding Greenwood’s formula (optional)

\[ \begin{aligned} \hat{\text{S}}(t) &=\prod_{t_i<t}x_i/Y_i\\ \text{S}(t)&=\prod_{t_i<t}p_i \end{aligned} \]

\[ \begin{aligned} \frac{\hat{\text{S}}(t)}{\text{S}(t)} &= \prod_{t_i<t} \frac{x_i}. {p_iY_i} \\ &= \prod_{t_i<t} \frac{\hat p_i}{p_i} \\ &= \prod_{t_i<t} {\left(1+\frac{\hat p_i-p_i}{p_i}\right)} \\ &\approx 1+\sum_{t_i<t} \frac{\hat p_i-p_i}{p_i} \end{aligned} \]

\[ \begin{aligned} \text{Var}\left(\frac{\hat{\text{S}}(t)}{\text{S}(t)}\right) &\approx \text{Var}\left(1+\sum_{t_i<t} \frac{\hat p_i-p_i}{p_i}\right) \\ &=\sum_{t_i<t} \frac{1}{p_i^2}\frac{p_i(1-p_i)}{Y_i} \\ &= \sum_{t_i<t} \frac{(1-p_i)}{p_iY_i} \\ &\approx\sum_{t_i<t} \frac{(1-x_i/Y_i)}{x_i} \\ &=\sum_{t_i<t} \frac{Y_i-x_i}{x_iY_i} \\ &=\sum_{t_i<t} \frac{d_i}{Y_i(Y_i-d_i)} \\ \therefore\text{Var}{\left(\hat{\text{S}}(t)\right)} &\approx \hat{\text{S}}(t)^2\sum_{t_i<t} \frac{d_i}{Y_i(Y_i-d_i)} \end{aligned} \]

Test for differences among the disease groups

Here we compute a chi-square test for assocation between disease group (group) and disease-free survival:

survdiff(surv ~ group, data = bmt)
#> Call:
#> survdiff(formula = surv ~ group, data = bmt)
#> 
#>                      N Observed Expected (O-E)^2/E (O-E)^2/V
#> group=ALL           38       24     21.9     0.211     0.289
#> group=Low Risk AML  54       25     40.0     5.604    11.012
#> group=High Risk AML 45       34     21.2     7.756    10.529
#> 
#>  Chisq= 13.8  on 2 degrees of freedom, p= 0.001

Cumulative Hazard

\[ \begin{aligned} {\lambda}(t) &\stackrel{\text{def}}{=}\text{p}(T=t|T\ge t)\\ &= \frac{\text{p}(T=t)}{\text{P}(T\ge t)}\\ &= -\frac{\partial}{\partial t}\text{log}{\left\{\text{S}(t)\right\}} \end{aligned} \]

The cumulative hazard (or integrated hazard) function is

\[{\Lambda}(t)\stackrel{\text{def}}{=}\int_0^t {\lambda}(t) dt\]

Since \({\lambda}(t) = -\frac{\partial}{\partial t}\text{log}{\left\{\text{S}(t)\right\}}\) as shown above, we have:

\[ {\Lambda}(t)=-\text{log}{\left\{\text{S}(t)\right\}} \]

So we can estimate \({\Lambda}(t)\) as:

\[ \begin{aligned} \hat {\Lambda}(t) &= -\text{log}{\left\{\hat{\text{S}}(t)\right\}}\\ &= -\text{log}{\left\{\prod_{t_i < t}{\left[1-\frac{d_i}{Y_i}\right]}\right\}}\\ &= -\sum_{t_i < t}\text{log}{\left\{1-\frac{d_i}{Y_i}\right\}} \end{aligned} \]

This is the Kaplan-Meier (product-limit) estimate of cumulative hazard.

Example: Cumulative Hazard Curves for Bone-Marrow Transplant (bmt) data

[R code]

autoplot(
  fun = "cumhaz",
  km_model1,
  conf.int = FALSE,
  ylab = "Cumulative hazard (disease-free survival)",
  xlab = "Time since transplant (days)"
) +
  theme_bw() +
  theme(legend.position = "bottom")

Since \(\text{S}(t)=\text{exp}\left\{-{\Lambda}(t)\right\}\), the Nelson-Aalen cumulative hazard estimate can be converted into an alternate estimate of the survival function:

\[ \begin{aligned} \hat S_{NA}(t) &= \text{exp}\left\{-\hat H_{NA}(t)\right\}\\ &= \text{exp}\left\{-\sum_{t_i < t}\frac{d_i}{Y_i}\right\}\\ &= \prod_{t_i < t}\text{exp}\left\{-\frac{d_i}{Y_i}\right\}\\ \end{aligned} \]

Compare these with the corresponding Kaplan-Meier estimates:

\[ \begin{aligned} \hat H_{KM}(t) &= -\sum_{t_i < t}\text{log}\left\{1-\frac{d_i}{Y_i}\right\}\\ \hat S_{KM}(t) &= \prod_{t_i < t}\left[1-\frac{d_i}{Y_i}\right] \end{aligned} \]